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SUMMARY 


A new inverse boundary layer method is developed and applied to incompressible 
flows with laminar separation and reattachment. Test cases for two-dimensional 
flows are computed and the results are compared with those of other inverse 
methods. One advantage of the present method is that the calculation of the 
inviscid velocities may be determined at each marching step without having to 
iterate. 

The inverse method was incorporated with the direct method to calculate the 
incompressible, conical flow over a slender delta wing at incidence. The location 
of the secondary separation line on the leeward surface of the wing is determined 
and compared with experiment for a unit aspect ratio wing at 20.5 deg incidence. 
The viscous flow in the separated region was calculated using prescribed 
skin-friction coefficients. 


INTRODUCTION 


The flow field over slender, highly swept delta wings at moderate incidence is 
dominated by the presence of large counter-rotating, leading-edge vortices as shown 
in fig. 1. As the flow moves spanwise towards the leading edge, the adverse 
pressure gradient caused by the leading-edge vortices causes the boundary layer to 
separate along a secondary separation line, indicated in fig. 2. 

A direct boundary layer method (one in which the external pressure is pre- 
scribed from an inviscid calculation or experiment) may be used to determine the 
location of the secondary separation line (ref. 1). However, in order to continue 
the solution from the secondary separation line to the leading-edge using boundary 
layer theory, an inverse method (one in which the wall shear or displacement thick- 
ness is specified) must be employed. 

Inverse methods have been used by numerous authors since the early work of 
Catherall and Mangier (ref. 2). Catherall and Mangier used a prescribed displace- 
ment thickness distribution to drive a boundary layer method in which the external 
pressure was determined as part of the solution. In this manner, they were able to 
obtain a regular solution at separation. However, their numerical scheme 
developed instabilities in the reversed flow region and the integration was con- 
tinued only by reducing the convergence criteria at each marching step. This 
problem of reversed-f low velocity profiles led directly to the FLARE approximation 
of Reyhner and Flugge-Lotz (ref. 3). In the FLARE approximation the streamwise 
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convection of momentum is set equal to zero in the reversed flow region and the 
calculations proceed with the usual forward marching procedure. 

For two-dimensional flows, Klineberg and Steger (ref. 4), Carter (ref. 5) and 
Cebeci et al. (ref. 6) have developed globally iterative schemes for flows with 
separation and reattachment. In those methods either the wall shear or displace- 
ment thickness is specified and the pressure is obtained iteratively using succes- 
sive-under-relaxation schemes. Semi-inverse methods have been developed for inter- 
acting an inverse boundary layer method with an inviscid solution by Le Balleur 
(ref. 7), Carter (ref. 8), Kwon and Pletcher (ref. 9) and Veldman (ref. 10) among 
others. In these methods, the inviscid calculation proceeds in the direct manner 
with the viscous calculation performed in the inverse mode. Thus, an edge velocity 
is determined in both calculations which must be the same after convergence. 

In three-dimensional flow calculations, several possible combinations of the 
viscous parameters could be used to drive an inverse method (e.g., CFX, CFY or DTX, 
DTY). However, Edwards and Carter (ref. 11) have shown that specifying the three- 
dimensional displacement surface and the component of vorticity normal to the 
surface leads to an elliptic set of equations and departure solutions for forward 
marching schemes. 

In this paper, a new inverse boundary layer method is developed for separated 
flows. The method is non-iterative based on a predictor-corrector linearization 
of the discretized governing equations. For two-dimensional flows, the new method 
is used to compute the test case of Klineberg and Steger (ref. 4) for a specified 
wall shear and that of Carter (ref. 5) for the displacement thickness prescribed. 
The method is also applied to the incompressible, conical, laminar boundary layer 
flow on the leeward surface of a slender delta wing at incidence. The inviscid 
solution for the delta wing is determined using the Free Vortex Sheet (FVS) code 
which was orginally developed at Boeing (ref. 12) and significantly enhanced by 
Luckring and others (refs. 13, 14) at NASA Langley. The viscous solution for the 
delta wing is determined in the direct mode until the secondary separation line is 
encountered (CFY < 0), then the calculations are continued in the inverse mode to 
the leading edge by specifying the skin-friction coefficients to be constant and 
equal to their values at separation (CFX > 0, CFY < 0). The calculations are per- 
formed on a unit aspect ratio wing at 20.5 deg incidence which corresponds to the 
experiment of Hummel (ref. 15). 


ANALYSIS 


Governing Equations 

For steady, incompressible flows, the three-dimensional laminar boundary layer 
equations in non-dimensional Cartesian coordinates are: 
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The velocity components u, v, and w ate in the x, y, and z-directions , respec- 
tively, where x is in the streamwise direction, y is spanwise and z is normal to 
the x-y plane. The Reynolds number has been removed from eqs. (l)-(3) by defining 

x = x/L, y = y/L, z “ /Re z/L, u = u/U os v = v/U 0 and w - » / Re w/U Q . The super- 
script bars indicate dimensional quantities and Re = U 0 L/v. The complete derivation 
of eqs. ( 1 )— ( 3) is provided in reference 16. The boundary conditions on eqs. 

( 1 ) — ( 3 ) are: 


at z = 0: u = v = w = 0 


and as z + °°: u*u, v+v 

e e 

Two-dimensional Flow 

For two-dimensional flow, v = 0 and eqs. (1) and (3) are transformed by 
defining, 


(4) 


1/2 


X = x 


= (— ) 

'■x ' 


along with a stream function, ip , where 
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Define a parameter f, such that 
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and (7) in eq. (6) yields, 
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With eqs. (5)-(8), eq. (2) may be written as, 


Using eqs. (5) 
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The boundary conditions for eq. (9) are 

p = 0 : f = F = 0 , n = q e ; F=1.0 


( 10 ) 
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Conical Flow Over Delta Wings 

For the delta wing a body-oriented coordinate system is defined through the 
transformation, 

X = x , Y = (y/x)tan A , Z = z/x (12) 

where A is the sweep angle of the delta wing. Define transformed velocity com- 
ponents by, 

U = u , V = v tan A-uY, W = w - uZ (13) 

Using eqs. (12) and (13) in eq. (1) it becomes, 

Iv (X 2 U) + It 7 (XV) + — (XW) = 0 (14) 


Next define two stream functins, ijj and 4> as 
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and a boundary layer variable, n> by 

n - (u x) 1/2 z 

e 

Introduce two functions f(X,Y,n) and g(X,Y,n) where, 
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For conical inviscid flow, fluid properties do not vary along rays; therefore, 
for any flow parameter, H, 3H/3X = 0. Then using eqs. (12), (16) and (17) in eq« 

(15), the following are obtained: 
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Equation (18) may be written as the vector equation, 
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With eqs. ( 12 )-( 20 ) , eqs. (2) and (3) may be combined and written as, 
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The boundary conditions for eqs. (19) and (21) are: 
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Equations (19) and (21) may be solved in the direct mode, using prescribed inviscid 
velocities, to calculate viscous properties such as skin-friction coefficients and 
integral thicknesses® However, in the inverse mode, the inviscid velocities are 
unknown functions and must therefore be eliminated from the left-hand side of eq. 
(21)® This is accomplished through the additional transformation, 
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Using eq. (25) in eq. (21), the latter becomes, 
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Numerical Procedure 


The governing equations are differenced using the half-implicit finite- 
difference scheme of Matsuno (ref. 17). The scheme is second-order accurate and 
unconditionally stable, and was demonstrated by Woodson and DeJarnette (ref. 18) to 
yield accurate numerical results when compared to the exact solution of the three- 
dimensional boundary layer equations for parabolic flow over a moving flat plate. 


For fully three-dimensional problems the scheme has the advantage that the 
crossflow derivative formulas are formed independent of the direction of the cross- 
flow. Further, the scheme is non-iterative based on a predictor-corrector linear- 
ization, a convenient feature for inverse boundary layer methods. For any flow 
parameter, H, the notation Hj^ = H(Y*,ri^) is used where, 


Y* = Y* + AY* , j = 1,2,... , JMAX 
j+1 j 

n n 

k+l = k + An , k = 1,2,...,KMAX 


(28) 


Central difference operators are defined by: 


A n H j,k = < H d > k+1 " Hj ,k-l)/An 


(29) 
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(30) 


A Y*H j+1/2>k - (Hj+1 ,k - 

5 n H j ,k-l/2 = (Hj,k “ H j ,k™ 1 )/ A n (31) 

6 n H j s k ~ ( H j ,k+l “ 2Hj jk + Hj >k _i)/Ar|2 (32) 

and a backward difference operator for the predictor stage by, 

VY*Hj+l/2,k “ 2(H j+l/2 , k ~ H jjk )/AY* (33) 

With the operators defined above, eqs. (19) and (26) become at the predictor stage, 
respectively, 

6 q h j+l/2,k-l/2 = (Hj+l/2 ,k + H j+l/2,k-l)/2 (34) 

^n^j+l/2,k + a j ,k A n H j ,k + ^Y*Sj+l/ 2 ,k A nHj jk - Gj j k VY*Hj+i/ 2 }k = lj,k (35) 

at the corrector stage they are, 

<$n h j+l,k-l/2 " ( H j+l,k + H j+l,k-l)/2 (36) 

2 

<Sr|(Hj+i > k + Hj jk )/2 + Oj+i/2 ,k A TiHj+i/2 >k 

(37) 

+ A Y*gj+l/2,k A n H j+l/2,k " G j+l/2,k A Y* H j+l/2,k = r j+l/2,k 
Equations (34)-(37) may be written in the block tridiagonal matrix form, 

h k = h k-l + (®k + H k-l) (38a) 

" H k _! + B k H k - H k+1 + a k h k = D k (38b) 

wnere a k and B k are 2x2 coefficient matrices and D k is a vector. Equations (38) 
are linear tridiagonal matrix equations and may be solved using a block form of the 
Davis modified tridiagonal alogrithm (ref. 19). 
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Inverse Method 


The term D k in eq. (38b) contains terms with both $ u and g v as coefficients. 

However, in the inverse mode these parameters are unknowns which must be determined 
from the solution. Taking advantage of the fact that eq. (38) is linear, one may 
write, 
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where the 

unknowns 8 U and 

8v are given by eq. (27). 

The boundary conditions 


associated with eq. (39) are: 


at n = 0 : hj 1 ^ = = 0 , i = 0,1,2 


and at n = h e : H KMAX = 1 > H KMAX = 0 > i = 1,2 


(40) 


Carrying out the operations of eqs. (29) -(33), one may equate like powers of f3 u and 
8 V to obtain 


hk £) = h k-i + f 11 ( H ^ ±) “ 4-1) » 1 “ 0,1,2 (41a) 

- + B k 4 i} - H^+1 + a k h k l) = d£ i} , i « 0,1,2 (41b) 


This system of equations may be solved to obtain h^ 1 and H k for i = 0,1,2. The 
parameters U e and V e are then determined from either the specified wall shear 
stresses or integral thicknesses. A streamwise integral thickness is defined by, 
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Using eqs. (16) and (18a) in eq. (42) and integrating gives. 
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DTX = 


(43a) 
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“ /5£T - n e - f(n e ) 


In a similar manner, 
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With eq. (39c), eq. (43) may be written as, 
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Equations (44a) and (44b) are solved to obtain B u and g v . Now, returning to 
eq. (27) and taking second-order differences about the point j+1/2, 
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Next, consider if the wall shear were perscribed. The surface skin-friction 
coefficients are defined by 
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Using eq. (5) in eq. (46) one obtainSj 
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(47b) 
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Using second-order accurate difference expressions for the derivatives in eq. (47) 
yields , 
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Substituting eq. (39b) into eq. (48) gives two equations for the two unknowns, g u 
and g v . An advantage of the present method over some earlier methods (refs. 4-6) 
is that the solution for f? u and 6 V is obtained without column iteration, as opposed 
to using an under-relaxation scheme. For regions of reversed flow the FLARE 
approximation is made at the predictor stage, 
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while at the corrector stage. 
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Equation (49) insures diagonal dominance of the Davis modified algorithm in the 
reversed flow region. Note that eq. (49) reduces the accuracy of the finite- 
difference method to first order. Carter (ref. 5) added an artificial "time-like” 
term to the term in eq. (38b) to ensure unconditional diagonal dominance in his 
solution algorithm. For the cases computed thus far, no instabilities were 
encountered in the reversed flow region when evaluating the convective derivatives 
according to eq. (49). 


RESULTS AND DISCUSSION 


Two-dimensional Flow 

Calculations using the inverse boundary layer method are first discussed for 
flow over a circular cylinder. The inviscid solution for flow over a circular 
cylinder is given from potential theory by, 

u e = 2sin 6 (50) 
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where 0 is the angle measured with respect to the x axis which defines the cylinder 
(0 < 0 < tt ) . The boundary layer was computed in the direct mode until separation 
was encountered at 0 = 104.5 deg and then the calculations were continued in the 
inverse mode by specifying the skin-friction coefficient to be constant and equal 

to its value one marching step beyond separation (CFX = 0 or slightly negative). 

Let f in eq. (8) and F in eq, (9) be written as. 
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+ mF k 

(51b) 


with corresponding boundary conditions, 


at 


n = 0 : fi l) = f{ 1) = 0 , i = 0,1 


and at n = b e * FKMAX = 1 . ^KMAX = 0 


(52) 


After substituting eq. (51) into eq. (9), two separate equations result; one for 
the variables with superscript (0) and another for superscript (1). After solving 
these two equations, the value of m can be calculated by substituting eq. (51) into 
eq. (48a). The values of m can then be used to numerically calculate u e at each 
marching station. The resulting edge velocity is shown in fig. 3. The edge 
velocity from the inverse calculation departs tangentially from the potential flow 
curve at the separation point and 6 = it is no longer a stagnation point. The 
inverse method could be used to calculate the edge velocity over the entire 
cylinder provided the wall shear or displacement thickness distribution is known. 

Klineberg and Steger (ref. 4) constructed a test case for a flow with separa- 
tion and reattachment by prescribing the skin-friction distribution as. 
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= T . 


X < 2, X > 6 


CFX = t [1 + oi(X-2) (X-6) ] , 2 < X < 6 


(53) 


where a is a given parameter. The resulting pressure gradient parameter, m, for 

a - 0.225 was calculated and is compared to that which was determined by Klineberg 
and Steger in fig. 4. Klineberg and Steger required between 400 to 800 iterations 
using a point successive-under-relaxation scheme for the calculation of m while the 
present results were obtained with one downstream pass . The corresponding displace- 
ment thickness and edge velocity for this case are given in fig. 5, and velocity 
profiles at three marching stations are shown in fig. 6. It was found that for 
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attached flow the profile for F^ was primarily positive, while for separated flow 

( 1 ) 

the majority of the profiles were negative. 

A second test case is that of Carter (ref. 5) in which the displacement thick- 
ness is perscribed according to the relation. 


DTX = ai + a 2 (X-X 0 ) + a 3 ( X-X Q ) 2 + a 4 (X-X 0 ) 3 , X < X Q 
DTX = a x + a 2 (X-X 1 ) + a 3 (X-X 1 ) 2 + a^X-X^ 3 , Xi < X < X 2 


(54) 


where X Q = 1.065, X 3 = 1.35, X 2 = 1.884 and the values of the a's are determined 
such that at X = X Q the value and slope of the displacement thickness match the 
Blasius flat plate distribution, and at X = X 3 , DTX reaches a maximum value. The 
case identified by Carter as Case B has a maximum displacement thickness of 8.6 and 
its distribution is presented in fig. 7a along with the resulting skin-friction 
distribution in fig. 7b. The calculated value of m for this case is compared with 
Carter's results in fig. 8a. Carter also developed an "approximate forward march- 
ing" technique in addition to his globally iterative method, which employs the 
FLARE approximaton. Both of his methods, however, used an under-relaxation scheme 
for the calculation of m. For the grid indicated in fig. 8, Carter's globally 
iterative scheme converged in 130 iterations and his forward marching procedure 
required an average of 41 column iterations at each marching step. The present 
method required only one downstream pass. The edge velocity is presented in fig. 
8b. 

This test case has a more extensive separated flow region than the Klineberg 
and Steger case with the approach to separation and reattachment much steeper. 

Note that m reaches a relative minimum just prior to separation and reattachment 
which satisfies Meksyn's criterion (ref. 20) for a regular solution. Velocity 
profiles at two x-locations are given in fig. 9. At the point X = 1.393, nearly 
half of the profile is in the reversed flow region, however, the magnitude of the 
reversed flow velocity is small compared to the edge value. For this test case, 
the maximum negative value of u is about one-tenth of u e . 


Conical Flow Over Delta Wings 


The inviscid solution for the Hummel delta wing at 20.5 deg incidence was 
determined using the Free Vortex Sheet (FVS) theory (refs. 12-14). The FVS code 
solves the linearized Prandtl-Glauret equation for potential flow. It represents 
the wing, wake, and rolled-up vortex sheets by continuous quadratic doublet sheet 
distributions and the vortex core by a line distribution of vortices. Detailed 
pressure distributions as well as overall forces and moments are predictd by the 
FVS code; however, no effort is made to model the secondary vortices in the 
inviscid calculation. 

The viscous calculations were begun in the direct mode beginning at the 
reattachment line (which was found to be located along the wing centerline for this 
case) and marched spanwise until the secondary separation line was encountered. 
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The secondary separation line is assumed to be the line along which CFY first goes 
through zero. The solution was then continued in the inverse mode by specifying 
the skin-friction coefficients to be constant and equal to their values at separa- 
tion (CFY < 0 but CFX >0). The spanwise distribution of the skin-friction coeffi- 
cients is presented in fig. 10a with corresponding integral thicknesses in 10b. The 
secondary separation line was calculated to be along Y = 0.70. Hummel (ref. 15) 
observed from oil flow studies a secondary separation line along Y = 0.67. A 
discrepancy between the calculated and observed secondary separation line was 
expected, since the inviscid solution ignores completely the influence of the 
secondary vortices. Reference 1 showed that the secondary separation line could be 
calculated accurately with the direct boundary layer method when experimental 
pressure data were used. The momentum integral thicknesses and shape factors are 
shown in fig. 11. Both shape factors increase sharply in the vicinity of the 
secondary separation line with the spanwise shape factor remaining nearly constant 
afterwards while the streamwise shape factor decreases continually toward the 
leading edge. Velocity profiles at three spanwise stations are presented in fig. 

12. Both profiles indicate the spanwise thickening of the boundary layer and the G 
profiles beyond the secondary separation line show a small region of reversed flow 
near the surface. The inviscid velocities determined from the FVS code and those 
calculated with the inverse method are presented in fig. 13. As expected, the 
gradients in the inviscid velocities calculated from the inverse boundary layer 
method are much less steep than those of the inviscid solution calculated neglect- 
ing the boundary layer. 

Research is in progress to interact the boundary layer solution with the 
inviscid solution. Reference 21 used a three-dimensional integral turbulent 
boundary layer method to solve the flow field over a delta wing at incidence. The 
viscous and inviscid solutions were coupled and the resulting pressure distribution 
showed some improvement over the inviscid results near secondary separation. 


CONCLUDING REMARKS 


A three-dimensional, direct boundary layer method was extended to the inverse 
mode for separated flows. It is a predictor-corrector finite-difference method in 
which the FLARE approximation is made to the streamwise momentum term in the 
preditor stage but a finite-difference quotient is used for the corrector stage. 

This method allows the Inviscid velocity to be calculated without iterations by 
marching into the separated flow region using prescribed skin-friction coefficients 
or integral thicknesses. 

Two-dimensional test cases for laminar separation with prescribed skin friction 
or displacement thickness were found to compare well with other methods. Solutions 
were also obtained using the direct and inverse modes for conical inviscid flow over 
a delta wing at incidence. The direct mode was used to the secondary separation 
line and then the inverse mode continued the solution into the separated flow region 
with both skin-friction coefficients maintained at their values at secondary 
separation. Further studies are needed for the inverse mode in fully three- 
dimensional flows and interacting the boundary layer solution with the inviscid 
solution. 
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SYMBOLS 


a k » ®k 

f X f Y* 

CFX.CFY 

Dk 

DTX.DTY 

f »g 

F,G 

H,h 

HX.HY 

JMAX 

KMAX 

L 

m 

MTX.MTY 

Re 

Re x 

u,v,w 

U,V,W 

U 0 

x,y»z 

X,Y,Z 

y* 


coefficient matrices in eq. (38) 

skin-f riction coefficients defined by eq. (46) 


C, /Re , C c */Re 
fx x’ fy* x 


vector in eq. (38) 
_* _* 



functions defined by eq. (17) 
velocity ratios given by eq. (18) 
vectors defined by eq. (20) 
MTX/DTX, MTY/DTY 


number of mesh points in spanwise direction 

number of mesh points across the boundary layer 

reference length, m 

parameter defined by eq. (10) 

momentum integral thicknesses 

freestream Reynolds number, U c L/v 

local Reynolds number, (U e X)Re 

non-dimensional velocity components in x, y, and z directions, 
respectively 

transformed velocity components given by eq. (13) 
freestream velocity, m/s 

non-dimensional Cartesian coordinates streamwise, spanwise, and normal 
to the wing surface, respectively 

transformed coordinates given by eq. (12) 

transformed spanwise coordinate given by eq. (25) 


128 



parameter in eq. (53) 


' u 
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-* 

6 X 


V dU 
e e 

2 dY 


U 


dV 


U e dY 

3 V — *53 u 

parameter defined by eq. (23) 

central difference operators given by eqs. (31), (32) 

0 f (i - 1 -)»; 


-* 

6y* 


A 

V 
T) 

e 

A 

V 

a 


0 I m (i 

v e 

central difference operators defined by eqs. (29), (30) 

backward difference operator defined by eq. (33) 

transformed normal coordinate defined by eqs. (5), (16) 

angle of rotation for circular cylinder 

sweep angle of the delta wing 

kinematic viscosity coefficient, m^/s 

parameter defined by eq. (22) 

stream functions given by eqs. (7), (15) 


Subscripts 

e edge value 

j ,k mesh point locations given by eq. ( 28 ) 

v viscous calculation 

w wall 


Barred parameters are dimensional quantities 
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Fig. 2. Upper surface flow geometry. 
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A 


Present method 
Klineberg & Steger 
(ref. A) 



Fig. 3. Inviscid velocity over a 
circular cylinder. 


(a) 


(b) 


Fig. 4. Skin-friction coefficient and pressure 
gradient parameter. 





(a) (b) (a) (b) 

Fig. 5. Displacement thickness and edge Fig. 6. Velocity profiles for Klineberg 
velocity. and Steger test case. 


— Present method > AX=0.026, Arr-0.1. 11 ”15.0 
& Carter 



(a) (b) (a) (b) 

Fig. 7. Displacement thickness and skin- Fig. 3. Pressure gradient parameter and 

friction coefficient. edge velocity. 
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(a) 


(b) 


Fig. 9. Velocity profiles for Carter test 
case. 


(a) 


(b) 


Fig. 10. Skin-friction coefficients and 
integral thicknesses. 



Fig. 11. Momentum integral thicknesses 
and shape factors. 


Fig. 12. Velocity profiles for the delta 
wing. 



Fig. 13. Inviscid velocity components. 
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